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Abstract 

The present article is devoted to the influence of sediment layers on the process of 
tsunami generation. The main scope here is to demonstrate and especially quan- 
tify the effect of sedimentation on vertical displacements of the seabed due to an 
underwater earthquake. The fault is modelled as a Volterra-type dislocation in an 
elastic half-space. The elastodynamics equations are integrated with a finite element 
method. A comparison between two cases is performed. The first one corresponds 
to the classical situation of an elastic homogeneous and isotropic half-space, which 
is traditionally used for the generation of tsunamis. The second test case takes into 
account the presence of a sediment layer separating the oceanic column from the 
hard rock. Some important differences are revealed. We conjecture that deforma- 
tions in the generation region may be amplified by sedimentary deposits, at least for 
some parameter values. The mechanism of amplification is studied through careful 
numerical simulations. 
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1 Introduction 



The primary application of this study is that of tsunami generation by the 
deformation of the sea bottom following an underwater earthquake. We do 
not explicitly compute the tsunami waves induced above the generation region. 
The coupling between solid and water motions was already performed in our 
previous work [DD09] and can be done again if necessary. Here we are mainly 
interested in the extreme amplitudes of the seabed displacements during the 
first minutes of a tsunamigenic earthquake. Recall that the free surface motion 
roughly follows these displacements. There are two fundamental reasons for 
this. The first one is that the rupture velocity of the seismic source, V, is 
much larger than the phase velocity of the tsunami, c. In practice, for seismic 
sources, V is of the order of 3 km/s, whereas c is typically less than 250 
m/s, even for the deepest ocean basins [OS03]. It means that the gravitational 
forces do not have enough time to change the shape of the free surface during 
the characteristic time of the seabed motion [BMR72]. The second reason is 
that water is assumed to be incompressible and shallow. Altogether it means 
that for our purpose we can restrict our attention to the motion of the ocean 
bottom. Profiles of the ocean free surface are not computed in this paper. 

The two fundamental reasons mentioned above are often used to justify the 
passive approach for tsunami generation where the static sea-bed displacement 
is simply translated to the free surface to generate the initial condition. Our 
previous investigations [DDK06, DD07, DD09] showed important differences 
between passive and active generations when the resulting wave is generated 
by a moving bottom. 

We would like to underline that the present study is a theoretical one at this 
stage and must be followed by a thorough study of the influence of the various 
parameters such as rupture velocity, focal depth, dip angle, slip angle, fault size 
(only one-dimensional - ID - faults are considered in this paper), sedimentary 
structure (the sedimentary layer is taken as a Poisson solid in this paper). 
Moreover, the far- field effect might be different from the near-field effect. We 
simply explore some aspects of the influence of the sedimentary layering on 
tsunami generation. We do not consider historical examples, even though it 
should be done in the future. Our goal is to present a framework for studying 
the process of tsunami generation. Recall that ten years ago, Synolakis et 
al. [SLCY97] were writing: "There is a lack of quantitative information on 
sediment layers overlying tsunamigenic faults and about how these layers affect 
directly the generation of tsunamis." Our study is a small step toward a better 
understanding of the role of sediments. 

The influence of sedimentary layering was already mentioned in some studies 
[Fuk79, Oka88, FWB06]. Let us comment on the various results obtained so 
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far. Both studies [Fuk79] and [Oka88] point out that fracturing through thick 
sediments produces large displacements in the source region but relatively 
small displacements in the far field. In the work by Okal [Oka88] , the influence 
of the sediment layer was studied in the framework of normal modes and 
interesting results were obtained for sources inside as well as outside the layer 
of sediments. In the present study we perform direct numerical simulations by 
solving the elastodynamics equations with a finite element method (FEM). 

To our knowledge, the most recent numerical study concerning the role of 
sedimentation in sub duction- zone thrust faults is [FWB06]. The scope of that 
paper was the long-term evolution of a typical subduction wedge. A quite so- 
phisticated thermo-mechanical modelling of the plate movement with realistic 
rheology was used. As in our study, the governing equations were solved with 
a two-dimensional FEM. The authors came to some important conclusions 
(see also [SB06] for interesting remarks on sediment layers). We would like to 
quote some of them since there is a connection with our results: 

"Our numerical simulations demonstrate that sedimentation stabilizes the 
underlying wedge, preventing internal deformation beneath the basin. Maxi- 
mum slip during great-thrust earthquakes tends to occur where sedimentary 
basins stabilize the overlaying wedge. The lack of deformation in these stable 
regions increases the likelihood of thermal pressurization of the subduction 
thrust, allows the fault to load faster, and allows greater healing of the fault 
between rupture events." 

In view of the above results, it is interesting to compare the distribution of 
sediment thickness in the world oceans (see Figure 1) with the seismic hazard 
map (see Figure 2). As one can see on Figure 1, the sediment thickness varies 
from to 20 km. Important deposits can be found along the Eastern coasts of 
America and the Western coasts of Africa. Fortunately there is no substantial 
seismic activity in these regions. But in the Bay of Bengal, and in particular 
in the Andaman sea, the situation is different. In this part of the world two 
factors are present simultaneously: an important seismic hazard and thick sed- 
iment layers. Unfortunately, we do not have reliable information on sediment 
thickness in the Mediterranean region. One can find different information in 
the literature. The estimates go from 25 m [HRR+04] to 1500 m or even more 
[EOS + 05]. So, it is difficult to draw any conclusions. In this study we try to 
understand what kind of implications this may have for tsunami generation 
processes. 

The present paper is organised as follows. In Section 2 we briefly describe 
the simplified mathematical model which represents the Earth crust. In the 
same section we also give some ideas about the discretisation procedure of 
the governing equations. Section 3 contains the description of two idealised 
test cases. Then, some results on the profile of the seafloor are presented. Fi- 



3 



Total Sediment Thickness of the World's Oceans & Marginal Seas 



ISO" 17Q' 160' 150' 140' 130" 120* 110' 100" 90' 




ISO" 170' 160' 150' 140' 130 J 120" 110' 100" 90' 80° 70' 60" 50' 40' 30' 20" 



Thickness in Meters 



10000 20000 



A digital lota] su clinic in tin tiim <Ui;il);isc for the world's 
Geophysical Data Oilier (NGDO .Marine Geology k 

(.-.-■! in ■ '■ < Diii.-iiiii Tt.t tl;H;i !>n(l(j,-(j ivilli :i s|i:li [iiu til' -j 

rainipiled from llirit imiH-jgili' .-Kmnvri: previously pi iblislied 
isopach limps; ocean drilling results, both ODP and DSDP; 
and seismic reQection profiles archived at NCDC as well as 
seismic data and i^ip.it 1 1 n i:l|,. -ii iNNc :is part or the IOCs 
Geological/Gi-oplivsii al Atlas ol tilt Pacific (liAPA) project. 



1) Age of the imderlyiiu; «nisL 

IJJTirrtiuiii- liiMirry ..("lilt- i>r-c:i[| rnisl 
.'! r .Sri in Ilifii] lli-|ld:i III 1 lireTIienl. 

it !N; LtiaiH. .mil Wiiiir.ii i-i itiiiit iu * .litres, iiml 
ri> I In: ii.ihiiv i 'I llii: sc-linii'iDinn | ti h i.' ssi's cLBnx'i 
sedimems lo depocenlcrs 



niay iwnat'Oull; ripiv-t iii lik bust- ■■! itii- sediments. These 

data are intended in pn»it|n ;i minimum inlue for the 

llii. l.i', ..- llii- ...t'ln i i. i.l in ii [iiirhimliir ijcogiaphnr regjon. 



http://www.ngdc.noaa.gov/mgg/sedthicrt/sedttiick.htnnl 



Figure 1. Total sediment thickness of the world's ocean and marginal seas (source: 
NOAA). 

nally, important conclusions and practical recommendations to tsunami wave 
modelers are given in Section 4. Some directions for future research are also 
outlined. 



2 Mathematical model and numerical method 



In this paper we use the same mathematical model as in our previous study 
[DD09]. Nevertheless, we give here a brief description of the model and refer to 
[Dut07, DD09, KDID07] for more details. We emphasize again that we consider 
a ID fault. In other words the fault has only one dimension, its length L. 

The fault is assumed to lie inside a linear elastic isotropic material. In the 
next section both homogeneous and inhomogeneous distributions of the Earth 
crust properties will be considered. 
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Figure 2. Global seismic hazard map (source: Swiss Seismological Service). 

We provide the 3D version of the governing equations. Let a represent the 
stress tensor. The displacement field u = (u,v,w)(x,y, z,t) satisfies the clas- 
sical elastodynamic equations issued from continuum mechanics [AR02]: 

V-£ = P^, (1) 

where p the material density. It is common in seismology to assume that the 
stress tensor is determined by Hooke's law through the strain tensor e = 
\ (Vu + V*u). Therefore 

g = A(V • u)L + 2/i£, (2) 

where A and p are the Lame coefficients. The coefficient p is the shear modulus. 
Thus, we come to the following linear elastodynamic problem: 

V ■ (A(V ■ u)L + ^(Vn + V l u)) = p-^. (3) 

As already pointed out, these coefficients can possibly depend on the spatial 
coordinates (x,y,z) (x,y: horizontal, z: vertical). The Lame coefficients can 
be expressed in terms of Poisson's ratio v and Young's modulus E as follows: 

2(iu Ev E 
X = -, 7T = — VT, n~T> A* 



2u (1 + - 2u)' 2(l + i/)' 

This remark immediately leads to the following property of solutions to the 
steady version of Equation (3): such solutions do not depend on Young's mod- 
ulus but only on Poisson's ratio. 
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The fault is modeled as a dislocation inside an elastic material. This type of 
model is widely used for the interpretation of seismic motion. A dislocation 
is considered as a surface (in three-dimensional problems) or a line (in two- 
dimensional problems) in a continuous medium where the displacement field 
is discontinuous. The displacement vector is increased by the amount of the 
Burgers vector b along any contour C enclosing the dislocation surface (or 
line), i.e. 

jdu = b. (4) 

G 

We let a dislocation run at speed V along a fault inclined at an angle S with 
respect to the horizontal. The rupture starts at the point x = and z = —d 
(it is supposed to be infinitely long in the transverse y— direction), propagates 
at constant rupture speed V for a finite time L/V in the direction 8 and stops 
at a distance L. Let ( be a coordinate along the dislocation line. On the fault 
located in the interval < £ < L, the slip is assumed to be constant. The rise 
time is assumed to be 0. We insist here that in real situations, the fault is two- 
dimensional and the rupture velocity usually refers to the velocity at which 
the fault ruptures in the transverse ?/— direction (see for example [Oka88]). 



2. 1 Discretization of the elastodynamics equations 



In order to apply the FEM we rewrite the governing equation (3) in the domain 
Q in variational form. One has 



J p^- -v dVL + J g(u) : Vv dVt= J f ■ v dS, W G V, 



d u 
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where V is the linear closed subspace of ^if 1 (fi)j and / is the loading applied 
to the Neumann boundary T^. This term is equal to zero in our computations, 
since the seabed is considered to be a free boundary in geophysics. 

In order to discretize the time derivative operator we apply a classical second 
order finite-difference scheme. We underline that the resulting method is fully 
implicit and has the advantage of being free of any CFL-type condition. In 
such problems implicit schemes become advantageous since the velocity of 
propagation of seismic waves is of the order of 3 — 4 km/s. After discretizing 
in time, one obtains the following variational form: 



J p- 2 t 2 v dn + J a{u^)) : w dn = J f-v ds, 



l(n+l) _ 2 jj(n) _|_ ifn-1) 

where the superscript denotes the time step number, e.g. = u(x, t n ). Then, 
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Figure 3. Typical mesh used in the numerical computations. The length scale is in 
meters. 

we apply the usual P2 finite-element discretization procedure. For the numer- 
ical computations, we used the freely available code FreeFem++ [HPH098]. 

Let us say a few words about the boundary conditions and the treatment of 
the dislocation in the program. As already stated, the seabed is assumed to 
be a free surface: 

q_ ■ n — f — 0, z = 0. 

The other boundaries are assumed to be fixed or in other words, we apply 
Dirichlet type boundary conditions u — 0. The authors are aware of the re- 
flective properties of this type of boundary conditions. But we take a compu- 
tational domain which is sufficiently large, so that the seismic waves do not 
reach the boundaries during the simulation time. This approach is not com- 
putationally expensive since we use adaptive mesh algorithms [HPH098] and 
in the regions far from the fault, element sizes are considerably bigger than in 
the fault vicinity. A typical mesh used in simulations is plotted on Figure 3. 

Now, let us discuss the implementation of the dislocation. Across the fault, 
the displacement field is discontinuous and satisfies the following relation: 

ip~(x, t) — u~(x, t) = b(x, t), (5) 

where the signs ± denote the upper and lower boundary of the dislocation sur- 
face, respectively. In order to satisfy the condition (5) we apply the following 
boundary conditions on the fault surface: 

b b 
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Figure 4. Test case with a homogeneous medium. 

Remark: Due to the presence of huge hydrostatic pressures in the crust, the 
two sides of the fault cannot detach physically. In any case this situation does 
not occur in nature. Mathematically it means that the Burgers vector b is 
tangent to the dislocation surface at each point. 



3 Numerical results 



The numerical method described in the previous section was validated in our 
previous study [DD09] by comparing numerical results with analytical results. 
In the present work we compare vertical displacements at the ocean bottom 
in two different situations. The first test case corresponds to the traditional 
modelling procedure where the Earth crust is assumed to be a homogeneous 
elastic material. It is schematically depicted on Figure 4. For example, the 
well known Okada solution 1 [Oka85], which is still widely used to construct 
initial conditions for various tsunami propagation codes, is based on these 
assumptions. 

In the second test case, we add a sediment layer of thickness h s on top of the 
previous configuration. This situation is depicted on Figure 5. Let us provide 
some comments on the model for sediments chosen in this study. In fact, what 
we call sediments here is an elastic layer which has the mechanical properties of 
sand according to [Mei94] . It means that porosity and the two-phase nature of 
this medium are neglected. These effects should be investigated in the future. 
As stated by Okal [Oka88] , detailed oceanic models of sedimentary structure 
have evidenced strong gradients of seismic velocities. The chosen sedimentary 

1 The original paper by Okada was published in 1985. In the Russian literature this 
solution was already known in 1978, after the publication of results by Gusiakov 
[Gus78]. Some particular cases of Okada solution were known even earlier [MS71, 
FB76]. 
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Figure 5. Test case with a sediment layer in between the oceanic column and the 
hard rock. 



parameter 


value 


Fault depth, d 


4000 m 


Dip angle, 5 


13° 


Fault length, L 


2000 m 


Slip along the fault, b 


10 m 


Fault propagation velocity, V 


2500 m/s 



Table 1 

Values of fault parameters used in this study. 



parameter 


sand 


granite 


Shear modulus fj,, Pa 


2 x 10 8 


30 x 10 9 


Poisson's ratio 


0.3 


0.27 


Shear wave velocity y/fi/p, m/s 


330 


3230 



Table 2 

Values of mechanical parameters for sand and granite. 

model is not extreme; much looser and weaker structures are encountered in 
nature. 

Let us now discuss the results. First we present static solutions corresponding 
to the test cases described above. Two solutions are plotted on Figure 6. In this 
case we take a sediment layer thickness h s equal to 600 m. In all figures of this 
section we plot the vertical displacement at the free surface of the Earth crust 
(or at the seabed, in other words). The values of the other parameters used in 
the computations are given in Tables 1 and 2. As the reader can see, there is no 
significant difference between the two solutions. In other words, sediments do 
not influence the static deformations due to a dislocation source. Physically, 
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Figure 6. Volterra dislocation source. Static solutions with (dashed line) and without 
(solid line) sediments. The thickness of the sediment layer is h s = 600 m. 

we can understand this situation, since in the static case the sand layer is just 
raised up by the deformed granite. Note that the parameters used here are not 
realistic. Indeed the length of the fault L is quite small and we repeat that the 
fault is ID. The slip b is extremely large but since at this stage the problem 
is linear, the absolute value of the slip is not that important. 

There is another rather mathematical explanation. In fact, as mentioned above, 
the steady Lame equations 2 do not depend on Young's modulus but only on 
Poisson's ratio. Not surprisingly the Okada solution has the same property 
since the analytical expressions contain the Lame coefficients in the combina- 
tion A/ (A + /i), which depends only on Poisson's ratio v. One can see in Table 
2 that Poisson's ratio v is almost the same for sand and granite. This is why 
the sediment layer does not have a strong effect on the steady solution. 

Now, let us consider the dynamical effects. The results are presented on Fig- 
ures 7-9. In these computations we take the same thickness of the sand layer 
as in the static case (h s = 600 m). It can be seen on Figure 7 that the deforma- 
tion in the homogeneous case is initiated earlier. This is to be expected since 
the shear wave velocity in the sand is almost ten times slower than in the gran- 
ite. Later, the deformation in the inhomogeneous case starts to evolve. It is 
surprising that it produces much bigger displacements - see Figure 8. In other 
words, taking into account the sediments increases considerably the seabed 
deformations. When time evolves, both solutions eventually reach comparable 

2 We assume that we neglect volume forces as well. So, the system of governing 
equations is homogeneous. 
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Figure 7. Dynamic sea-bed displacements at the beginning of the rupture process. 
The thickness of the sediment layer is h s = 600 m. 
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Figure 8. Dynamic sea-bed displacements. The solution, which takes into account 
the presence of the sediments, produces much bigger vertical displacements. The 
thickness of the sediment layer is h s = 600 m. 

amplitudes (see Figure 9). 

We performed other computations where the sediment layer thickness was 
reduced to 150 m. Results are presented on Figures 10-11. In this case, both 
solutions (homogeneous and inhomogeneous) evolve together and are almost 
indistinguishable within graphical accuracy. These results suggest to study 
the dependence of the vertical displacement amplitude on the sediment layer 
thickness. 



3.1 Sediment amplification factor 



In order to quantify the influence of sediments on the vertical seabed displace- 
ments, we introduce a new quantity S a that is called the sediment amplification 
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Figure 9. Dynamic sea-bed displacements at the end of the simulation. The thickness 
of the sediment layer is h s = 600 m. 
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Figure 10. Dynamic sea-bed displacements at the beginning. The thickness of the 
sediment layer is h s = 150 m. 
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(a) t = 3.0s 



(b) t = 3.5s 



Figure 11. Dynamic sea-bed displacements. End of the process. The thickness of the 
sediment layer is h s = 150 m. 
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factor. We give first the formula for S a and then explain our definition. 

Definition 1 Let us denote by Vq(x, t) and v s (x,t) the vertical displacements 
at the free surface in a homogeneous half-space and at the top of the sedi- 
ment layer respectively 3 . Then, the sediment amplification factor is defined 
as follows: 

max \ v s (x, t)\ 

£ _ PM) J 

a max\vo(x,t)\ 

(x,t) 

Let us provide some explanations. First of all, it is clear that we compare 
the values of two extreme amplitudes. The maximum is taken in both space 
and time, since both processes are not synchronised in time 4 . Finally, we 
substract one because we want the amplification factor to be equal to zero 
when sediments are absent. 

Once this quantity S a is defined, a parametric study must be performed. Here 
we are mainly interested in the dependence on the sediment layer thickness. 
But it is better to choose a dimensionless quantity. In this problem there 
are three lengths: the fault length L, the fault depth d and the sediment layer 
thickness h s . It is natural to choose the ratio h s /d as dimensionless parameter. 

We performed a lot of computations for different values of h s /d and obtained 
the curve shown on Figure 12. It leads to several comments. On the left, the 
curve starts from zero and it is expected since the sediment layer disappears at 
this extremity. Thus, its amplification is equal to zero as well. It is interesting 
that the amplification factor has a maximum in the vicinity of h s /d = 0.12. It 
means that there exists an optimal configuration when the sediment layer has 
its strongest effect. When we gradually increase the dimensionless parameter 
h s /d past the maximum, the amplification decreases. In the limit, one has to be 
careful as h s /d — > oo. Indeed one approaches an elastic half-space completely 
filled with sediment and dislocation theory may no longer be appropriate, 
especially if the material is loose. 



4 Conclusions and perspectives 

In the present paper we investigated the influence of sedimentary layering on 
displacements due to an earthquake. We showed that there is practically no 
effect in the case of static deformation. This is to be expected in the framework 

3 In the idealized situation of our test cases, it means that we evaluate the vertical 
displacements at z = 0. 

We saw on Figures 7-9 that the homogeneous solution evolves faster since sedi- 
ments slow down the properties of elastic wave propagation. 
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Figure 12. Dependence of the sediment amplification factor S a on the sand layer 
relative depth h s /d. 

of our model. Both curves can be superimposed up to graphical accuracy. On 
the other hand, dynamics makes a big difference. Our computations show that 
the vertical displacement amplitude can be amplified by a factor up to 1.7. We 
point out that there exists some kind of "optimal" sediment layer thickness, 
which provides the biggest amplification factor. Of course, this optimal value 
depends on various mechanical parameters. It can be estimated in each specific 
situation by similar numerical techniques. 

There is another predictible effect of the sediment layer. It slows down con- 
siderably the velocity of propagation of elastic and Rayleigh waves. In our 
simulations it is reflected by the fact that the maximum amplitude is reached 
much later than in the homogeneous case. 

We introduced a new quantity S a , that we called sediment amplification factor. 
This dimensionless quantity measures the relative increase of the vertical dis- 
placement amplitude with respect to the homogeneous half-space solution. The 
dependence of this quantity S a on the dimensionless thickness of the sediment 
layer was studied. We showed that there exists an optimal ratio h s /d « 0.125 
between sediment thickness and depth of the event which provides the biggest 
amplification factor S a ~ 0.7. 

It is of interest to see how these results apply to the 2004 Indian Ocean 
Earthquake, which produced a megatsunami with local runup greater than 
30 m at some locations [SK06]. According to Figure 1, we can estimate the 
sediment thickness h s in the generation region to be about 3 km. In [LKA + 05, 
NSS + 05], the centroid depth d was set to 25 km for all fault subdivisions. 
If we compute the ratio of these parameters, we obtain h s /d = 0.12. This 
value approximatively corresponds to the value (see Figure 12) which provides 



14 



the maximal sediment amplification factor. The natural question is: is it a 
coincidence? It is not possible to answer at this stage. Indeed the values we 
have used for the fault size are not realistic. Most likely there can be an 
amplification in real situations but not as large. 

The overall conclusion of this study is that one may have to revise the initial 
conditions used in some tsunami simulations. More precisely one has to take 
care of situations where the generation region contains sediment deposits. 
Most likely it was the case of the Boxing Day Tsunami of 2004 [SB06] . Several 
researchers had to take unphysically large values of the slip along the fault 5 
in order to generate a significant tsunami wave (see for example [IAK + 07]). If 
one takes sediments into account, this value can be reduced while producing 
the same wave amplitude. 

We finally outline some directions for future research in this field. First of all, 
the application of these techniques to real world events requires, of course, 
3D computations, even if we do not think that it will change qualitatively 
our results. On the other hand, the fracturing through the sediments should 
be further investigated since it was conjectured to provide a much bigger 
amplification factor [Fuk79, Oka88]. At the same time, the question of the 
influence of sediment porosity has not been addressed in the present study 
and is left for future investigations. We have the feeling that porosity may 
enhance sea-bed deformations in the near field but it should be checked by 
thorough computations. 

Finally it should be emphasized that the present study has focused on the 
deformation of the sea-bed. The question of the amplification of tsunami waves 
is a different one. Even if there is an effect in the near-field, it is not obvious 
that there is an effect in the far-field. The normal mode theory used by Okal 
[Oka88] provides an ideal framework to study both near-field and far-field 
effects. Moreover it is relatively easy to vary as many parameters as possible 
(focal depth, rupture velocity, dip angle, slip angle, thickness of sediment layer, 
fault size, sedimentary structure, etc). 
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